FreeFem++, a tool to solve PDEs numerically 

Georges Sadaka 



(N 



< 



LAMFA CNRS UMR 7352 
Univcrsitc de Picardie Jules Verne 
33, rue Saint-Leu, 80039 Amiens, France 
http : / /lamf a . u-picardie . f r/ sadaka/ 



o 

CN ^ georges.sadaka@u-picardie.fr 

March 3, 2013 



Abstract 



FreeFem++ is an open source platform to solve partial differential equations numerically, 
based on finite element methods. It was developed at the Laboratoire Jacques-Louis Lions, 
Universite Pierre et Marie Curie, Paris by Frederic Hecht in collaboration with Olivier 
Pironneau, Jacques Morice, Antoine Le Hyaric and Kohji Ohtsuka. 

The FreeFem++ platform has been developed to facilitate teaching and basic research 
t-h through prototyping. FreeFem++ has an advanced automatic mesh generator, capable 

of a posteriori mesh adaptation; it has a general purpose elliptic solver interfaced with 
fast algorithms such as the multi-frontal method UMFPACK, SuperLU . Hyperbolic and 
parabolic problems are solved by iterative algorithms prescribed by the user with the 
high level language of FreeFem++. It has several triangular finite elements, including 
discontinuous elements. For the moment this platform is restricted to the numerical 
^vq simulations of problems which admit a variational formulation. 



We will give in the sequel an introduction to FreeFem++ which include the basic of this 
software. You may find more information throw this link http: //www. f reef em. org/ff++. 
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1 Introduction 



FreeFem++ is a Free software to solve PDE using the Finite element method and it run on 
Mac, Unix and Window architecture. 

In FreeFem++, it's used a user language to set and control the problem. This language allows 
for a quick specification of linear PDE's, with the variational formulation of a linear steady 
state problem and the user can write they own script to solve non linear problem and time 
depend problem. 

It's a interesting tool for the problem of average size. It's also a help for the modeling in 
the sense where it allows to obtain quickly numerical results which is useful for modifying a 
physical model, to clear the avenues of Mathematical analysis investigation, etc ... 

A documentation of FreeFem++ is accessible on www.freefem.org/ff++, on the following 
link www. f reef em. org/f f ++/ftp/FreeFem++doc .pdf , you may also have a documentation in 
Spanish on the following link http://www.freefem.org/ff++/ftp/freefem++Spanish.pdf 

You can also download an integrated environment called FreeFem++-cs, written by Antoine 
Le Hyaric on the following link www.ann.jussieu.fr/~lehyaric/ffcs/install.php 

2 Characteristics of FreeFemH — |- 

Many of FreeFem++ characteristics are cited in the full documentation of FreeFem++, we 
cite here some of them : 

- Multi- variables, multi-equations, bi-dimensional and three-dimensional static or time 
dependent, linear or nonlinear coupled systems; however the user is required to describe 
the iterative procedures which reduce the problem to a set of linear problems. 

- Easy geometric input by analytic description of boundaries by pieces, with specification 
by the user of the intersection of boundaries. 

- Automatic mesh generator, based on the Delaunay-Voronoi algorithm [LucPir98]. 

- load and save Mesh, solution. 

- Problem description (real or complex valued) by their variational formulations, the write 
of the variational formulation is too close for that written on a paper. 

- Metric-based anisotropic mesh adaptation. 

- A large variety of triangular finite elements : linear, quadratic Lagrangian elements and 
more, discontinuous PI and Raviart-Thomas elements, ... 

- Automatic Building of Mass/Rigid Matrices and second member. 

- Automatic interpolation of data from a mesh to an other one, so a finite element function 
is view as a function of (x; y) or as an array. 

- LU, Cholesky, Crout, CG, GMRES, UMFPack sparse linear solver. 

- Tools to define discontinuous Galerkin finite element formulations PO, Pldc, P2dc and 
keywords: jump, mean, intalledges. 

- Wide range of examples : Navier- Stokes, elasticity, fluid structure, eigenvalue problem, 
Schwarz' domain decomposition algorithm, residual error indicator, ... 

- Link with other software : modulef, emc2, medit, gnuplot, ... 
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- Generates Graphic/Text /File outputs. 

- A parallel version using mpi. 



3 How to start? 

All this information here are detailed in the FreeFem++ documentation. 

3.1 Install 

First open the following web page 

http : / /www . f reef em . org/f f ++/ 

Choose your platform: Linux, Windows, MacOS X, or go to the end of the page to get the 
full list of downloads and then install by double click on the appropriate file. 

3.2 Text editor 

1. For Windows : 

Install notepadH — |- which is available at http://notepad-plus.sourceforge.net/uk/ 
site .htm 

- Open Notepad++ and Enter F5 

- In the new window enter the command FreeFem++ " $ (FULL_CURRENT_PATH) " 

- Click on Save, and enter FreeFem++ in the box "Name", now choose the short cut key 
to launch directly FreeFem++ (for example alt+shif t+R) 

- To add Color Syntax Compatible with FreeFem++ In Notepad++, 

- In Menu "Parameters "->"Conf iguration of the Color Syntax" proceed as fol- 
lows: 

- In the list "Language" select C++ 

- Add "edp" in the field "add ext" 

- Select "INSTRUCTION WORD" in the list "Description" and in the field "supple 
mentary key word", cut and past the following list: 

PO PI P2 P3 P4 Pldc P2dc P3dc P4dc RTO RT1 RT2 RT3 RT4 macro plot intld 
int2d solve movemesh adaptmesh trunc checkmovemesh on func buildmesh square 
Eigenvalue min max imag exec LinearCG NLCG Newton BFGS LinearGMRES catch 
try intalledges jump average mean load savemesh convect abs sin cos tan atan asin 
acos cotan sinh cosh tanh cotanh atanh asinh acosh pow exp log loglO sqrt dx dy 
endl cout 

- Select "TYPE WORD" in the list "Description" and ... " "supplementary key word", 
cut and past the following list 

mesh real fespace varf matrix problem string border complex ifstream ofstream 

- Click on Save & Close. Now nodepad++ is configured. 

2. For MacOS : 

Install Smultron which is available at http://smultron.sourceforge.net. It comes 
ready with color syntax for .edp file. To teach it to launch FreeFem++ files, do a "command 
B" (i.e. the menu Tools/Handle Command/new command) and create a command which 
does 



4 



/usr/local/bin/FreeFem++-CoCoa 7 °/„p 



3. For Linux : 

Install Kate which is available at ftp://ftp.kde.Org/pub/kde/stable/3.5.10/src/ 
kdebase-3 . 5.10. tar . bz2 

To personalize with color syntax for .edp file, it suffices to take those given by Kate for 
C++ and to add the keywords of FreeFem++. Then, download edp.xml and save it in the 
directory " .kde/ share/apps/katepart /syntax" . 

We may find other description for other text editor in the full documentation of FreeFem++. 

3.3 Save and run 

All FreeFem++ code must be saved with file extension .edp and to run them you may double 
click on the file on MacOS or Windows otherwise we note that this can also be done in terminal 
mode by : FreeFem++ mycode . edp 



4 Syntax and some operators 
4.1 Data types 

In essence FreeFem++ is a compiler: its language is typed, polymorphic, with exception and 
reentrant. Every variable must be declared of a certain type, in a declarative statement; each 
statement are separated from the next by a semicolon " ; " . 

Another trick is to comment in and out by using the "//" as in C++. We note that, we can 
also comment a paragraph by using "/ * paragraph */" and in order to make a break during 
the computation, we can use "exit(O) ;". 

The variable verbosity changes the level of internal printing (0, nothing (unless there are 
syntax errors), 1 few, 10 lots, etc. ...), the default value is 2 and the variable clockO gives the 
computer clock. 

The language allows the manipulation of basic types : 



current coordinates : x, y and z; 

. . d d d d d 

current differentials operators : dx= — — , dy= — , dz= — dxy= — — , dxz= — — . 

ox oy oz oxy oxz' 

d d d d 

dyz= - — , dxx= - — , dyy= - — and dzz= — — ; 

oyz oxx oyy ozz 

integers, example : int a=l; 

reals, example : real b=l. ; (don't forget to put a point after the integer number) 
complex, example : complex c=l.+3i; 
strings, example : string test="toto"; 

arrays with real component, example: real [int] V(n) ; where n is the size of V, 
arrays with complex component, example: complex [int] V(n) ; 
matrix with real component, example: real [int , int] A(m,n) ; 
matrix with complex component, example: complex [int , int] C(m,n) ; 
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- bidimensional (2D) finite element meshes, example : mesh Th; 

- 2D finite element spaces, example : fespace Vh(Th,Pl) ; // where Vh is the Id space 

- threedimensional (3D) finite element meshes, example : mesh3 Th3; 

- 3D finite element spaces, example : fespace Vh3(Th3,P13d) ; 

- intld(Th,D ( u*v ) = J u ■ v dx where r C R; 

- int2d(Th) ( u*v ) = / u- v dxdy where Q C 1R 2 ; 

Jn 

- int3d(Th) ( u*v ) = / u- v dxdydz where C M. 3 . 

Jn 

4.2 Some operators 

We cite here some of the operator that are defined in FreeFem++: 

+, -, *, /, ~, 
<, >, <=, >=, 

k, I , // where a & b = a and b, a I b = a or b 



4.3 Manipulation of functions 

We can define a function as : 

- an analytical function, example : func uO=exp(-x A 2-y A 2) ,ul=l . *(x>=-2 & x<=2); 

- a finite element function or array, example : Vh uO=exp(-x A 2-y A 2) ;. 

We note that, in this case uO is a finite element, thus uO [] return the values of uO at each 
degree of freedom and to have access to the i element of uO [] we may use uO [] [i] . 
We can also have an access to the value of uO at the point (a,b) by using uO(a,b); 

- a complex value of finite element function or array, example : Vh<complex> uO=x+li*y; 

- a formal line function, example : func real g(int a, real b) { ; return a+b;} 

and to call this function for example we can use g(l,2). 

We can also put an array inside this function as : 

func real f(int a, real [int] U){ 
Vh NU; 
NU []=U; 
return a*NU; 

} 

Vh U = x , FNU = f (5 ,U[] ) ; 

- macro function, example : macro F(t ,u, v) [t*dx(u) ,t*dy (v)] //, notice that every 
macro must end by "//", it simply replaces the name F(t,u,v) by [t*dx(u) ,t*dy(v)] 
and to have access only to the first element of F, we can use F(t ,u, v) [0] . 

In fact, we note that the best way to define a function is to use macro function since in 
this example t,u and v could be integer, real, complex, array or finite element, ... 
For example, here is the most used function defined by a macro : 

macro Grad (u) [dx (u) , dy (u) ] // in 2D 
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macro Grad (u) [dx (u) , dy (u) , dz (u) ] // in 3D 

macro div (u , v) [dx (u) +dy ( v) ] // in 2D 

macro div (u , v , w) [dx (u) +dy (v) +dz (w) ] // in 3D 



4.4 Manipulation of arrays and matrices 

Like in matlab, we can define an array such as : real [int] U=l : 2 : 10 ; which is an array of 
5 values U[i]=2*i+1; i=0 to 4 and to have access to the element of U we may use U(i). 

Also we can define a matrix such as real [int , int] A=[ [1,2,3] , [2,3,4] ]; which is 
a matrix of size 2x3 and to have access to the (i, j)^ element of A we may use A(i, j). 

We will give here some of manipulation of array and matrix that we can do with FreeFem++: 

real [int] ul = [1 , 2 , 3] , u2 = 2 : 4 ; // defining ul and u2 

real ulpu2 = ul ' *u2 ; // give the scalar product of ul and u2 , here 

^►ul ' is the transpose of ul ; 
real [int] uldu2 = ul . / u2 ; // divided term by term 
real [int] ulmu2 = ul . *u2 ; // multiplied term by term 
matrix A=ul*u2 J ; // product of ul and the transpose of u2 
matrix <complex > C=[ [1 , 1 i] , [l + 2i , . 5*1 i] ]; 
real trA = trace ( [1 , 2 , 3] * [2 , 3 , 4] ' ) ; // trace of the matrix 
real detA = det ( [ [1,2], [-2,1] ]); // just for matrix lxl and 2x2 



4.5 Loops and conditions 

The for and while loops are implemented in FreeFem++ together with break and continue 
keywords. 

In for-loop, there are three parameters; the INITIALIZATION of a control variable, the 
CONDITION to continue, the CHANGE of the control variable. While CONDITION is true, 
for-loop continue. 

for (INITIALIZATION; CONDITION; CHANGE) 
{ BLOCK of calculations } 

An example below shows a sum from 1 to 10 with result is in sum, 
int sum=0; 

for (int i=l; i<=10; i++) 
sum += i ; 

The while-loop 

while (CONDITION) { 

BLOCK of calculations or change of control variables 

} 
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is executed repeatedly until CONDITION become false. The sum from 1 to 5 can also be 
computed by while, in this example, we want to show how we can exit from a loop in midstream 
by break and how the continue statement will pass the part from continue to the end of the 
loop : 

int i = 1 , sum = ; 
while (i<=10) { 

sum += i ; i ++ ; 

if (sum>0) continue; 

if (i==5) break ; 

} 



4.6 Input and output data 

The syntax of input/output statements is similar to C++ syntax. It uses cout, cin, endl, 
« and » : 

int i ; 

cout << " std-out" << endl; 
cout << " enter i= ? " ; 
cin >> i ; 
Vh uh=x+y; 

ofstream f ( " toto . txt " ) ; f << uh [] ; // to save the solution 
ifstream f (" toto . txt ") ; f >> uh [] ; // to read the solution 

We will present in the sequel, some useful script to use the FreeFem++ data with other 
software such as ffglut, Gnuplot 1 , Medit 2 , Matlab 3 , Mathematical Visit 5 when we save 
data with extension as .eps, .gnu, .gp, .mesh, .sol, .bb, .txt and .vtu. 

For ffglut which is the visualization tools through a pipe of FreeFem++, we can plot the solution 
and save it with a .eps format such as : 

plot (uh , cmm="t = "+t + " ; I I u I I _L ~ 2= " + N0RML2 [kk] , f ill = true , value = 
*true , dim=2) ; 

For Gnuplot, we can save the data with extension .gnu or .gp such as : 
{ ofstream gnu (" plot . gnu ") ; // or plot . gp 

//ofstream gnu (" plot ." + 1000 + k ". gnu ") ; // to save the data 
for (int i =0 ; i <=n ; i ++) 

gnu<<xx [i] << " " <<yy [i] <<endl ; // to plot yy[i] vs xx[i] 

} 

exec("echo 'plot \"plot.gnu\" w lp \ 
pause 5 \ 

quit ' | gnuplot " ) ; 



1. http://www.gnuplot.info/ 

2. http : //www . ann . jussieu . f r/~f rey/ software . html 

3. http : //www.mathworks . f r/products/matlab/ 

4. http : //www. wolf ram. com/mathematica/ 

5. https : //wci . llnl .gov/ codes/visit/ 
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teO ;||u||_L*2= 1.25 193 




■0.138515 

■0.19393 

■0.2*326 

■0.304732 

■0.360138 

■0.415544 

■0.47095 

■0.526356 

■0.581761 

■0.637167 

■0.692573 

■0.747979 

■0.803385 

■0.858791 

■0.914197 

■0.969602 

■1.02501 



Figure 1: Visualising of the solution using ffglut 



For Medit, we can save the data with extension .mesh and .sol such as 

load "medit " 
int k=0; 

savemesh (Th , "solution . "+(1000+k)+" .mesh") ; 

save sol ( " solution . " + (1000+k) +" . sol " , Th , uh) ; 

medit (" solution ", Th , uh) ; // to plot the solution here 

k + =l ; 



And then throw a terminal, in order to visualize the movie of the first 11 saved data, we can 
type this line : 

ffmedit -a 1000 1010 solution . 1000 . sol 

Don't forget in the window of Medit, to click on "m" to visualize the solution! 
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3-:03E-0: 





Figure 2: Visualising of the solution using Medit 

For Matlab, we can save the data with extension .bb such as : 

{ ofstream f ile ( " solut ion . bb " ) ; 

file << "2 1 1 "<< Vh.ndof << " 2 \n"; 
for (int j =0 ; j <Vh . ndof ; 

file << uh[][j] << endl ; 

} 



And in order to visualize with Matlab, you can see the script made by Julien Dambrine at 
http : / /www . downloadplex . com/Publishers/ Julien-Dambrine/Page- 1-0-0-0-0 . html. 




Figure 3: Visualising of the solution and the mesh using matlab 
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For Mathematica, we can save the data with extension .txt such as : 
int k=0; 

{ ofstream f f ( " uhsol . " + ( 1000 + k) + 11 . txt " ) ; 
for (int i=0 ; i<Th . nt ; i++) { 

for (int j=0; j <3; 

ff <<Th [i] [j] . x<<" "<< Th[i][j].y«" " <<uh [] [ Vh ( i , j ) ] << endl ; 

ff <<Th [i] [0] . x<<" "<< Th [i] [0] . y<<" " <<uh [] [Vh (i , 0) ] <<" \n" ; 

} 
} 

k+=l ; 



Aullior : i » i u i ■■ i . S ADA KA — LAMFA — Amiens 
BBM-BBM Motilon on a HARBOR M 1 = 0» |aec] 




Figure 4: Visualising of the solution using Mathematica 



For Visit, we can save the data with extension .vtu such as : 

load "iovtk" 
int k=0; 

int[int] f f order2= [1 , 1 , 1] ; 

savevtk ("solution. "+(1000+k)+" . vtu" , Th , uhl , uh2 , order =fforder2 , 

**dataname="UHl UH2 " , bin=true); 
k + =l ; 
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Figure 5: Visualising of the solution using visit 



5 Construction of the domain Q 

We note that in FreeFem++ the domain is assumed to described by its boundary that is on 
the left side of the boundary which is implicitly oriented by the parametrization. 
Let O be the rectangle defined by its frontier dfl = [—5,5] x [—1, 1] where his vertices are 
A(-5, -1), 5(5, -1), C(5, 1) and D(-5, 1), so we must define the border AB, BC, CD and DA 
of dQ by using the keyword border then the triangulation Th of Q is automatically generated 
by using the keyword buildmesh. 

real Dx=.2; // discretization space parameter 
int aa=-5 , bb=5 , cc=-l , dd=l ; 



border 


AB 


(t 


= aa , 


bb) {x 


= t ;y = 


cc 


label 


= i;} 


border 


BC 


(t 


= cc , 


dd) {x 


= bb;y = 


t 


label 


= 2;} 


border 


CD 


(t 


= bb , 


aa) {x 


= t ;y = 


dd 


label 


= 3;} 


border 


DA 


(t 


= dd, 


cc ) {x 


= aa;y = 


t 


label 


= 4;} 



mesh Th = buildmesh( AB ( floor ( abs (bb-aa) /Dx) ) + BC ( f loor ( abs (dd- 
k»cc)/Dx)) + CD (floor (abs (bb-aa) /Dx) ) + DA ( f loor ( abs ( dd- cc ) /Dx ) 

*0 ); 

plot( AB (floor (abs (bb-aa) /Dx) ) + BC ( floor ( abs (dd-cc ) /Dx) ) + CD ( 
**f loor (abs (bb-aa) /Dx) ) + DA (floor (abs (dd-cc) /Dx) ) ); //to see 
the border 

plot ( Th , ps = " mesh . eps " ) ; // to see and save the mesh 

The keyword label can be added to define a group of boundaries for later use (Boundary 
Conditions for instance). Boundaries can be referred to either by name ( AB for example) or 
by label ( 1 here). 
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A B 

Figure 6: Plot of the border (left) and the mesh (right) 

Another way to construct a rectangle domain with isotropic triangle is to use : 

mesh Th = square (m , n , [x , y] ) ; // build a square with m point on x 

^direction and n point on y direction 
mesh Thl =movemesh (Th , [x+1 , y *2] ) ; // translate the square 

>■*•] , 1 [*] , 1 [ to a rectangle ]1,2[*]0,2[ 
savemesh (Thl ," Name . msh ") ; // to save the mesh 
mesh Th2 ( 11 mesh . msh ") ; // to load the mesh 



label=3 




label=1 



Figure 7: Boundary labels of the mesh by square (10, 10) 

We can also construct our domain defined by a parametric coordinate as: 

border C(t=0,2*pi){ x= cos (t ); y= sin (t ); label =1} 
mesh Mesh_Name=buildmesh (C (50) ) ; 



13 



Figure 8: mesh Th by build (C (50) ) 



To create a domain with a hole we can proceed as: 

border a(t=0,2*pi){ x=cos(t); y= sin (t ); label = 1 ; } 

border b(t=0,2*pi){ x = . 3 + . 3* cos (t ) ; y = . 3* sin (t ) ; label =2 ; } 

mesh Thwithouthole= buildmesh (a (50) +b (+30) ) ; 

mesh Thwithhole = buildmesh (a (50) +b ( -30) ) ; 

plot ( Thwithouthole ,wait=l ,ps="Thwithouthole . eps") ; 

plot (Thwithhole , wait=l ,ps=" Thwithhole . eps") ; 




Figure 9: mesh without hole Figure 10: mesh with hole 



6 Finite Element Space 

A finite element space (F.E.S) is, usually, a space of polynomial functions on elements of Th, 
triangles here, with certain matching properties at edges, vertices, ... ; it's defined as : 

f espace Vh ( Th , PI ) ; 
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As of today, the known types of F.E.S. are: PO, P03d, PI, P13d, Pldc, Plb, Plb3d, P2, 
P23d, P2b, P2dc, P3, P3dc, P4, P4dc, Morley, P2BR, RTO, RT03d, RTOOrtho, 
Edge03d, Pine, RT1, RTlOrtho, BDM1, BDMlOrtho, TDNNS1; where for example: 

P0,P03d piecewise constant discontinuous finite element (2d, 3d), the degrees of freedom are 
the barycenter element value. 

P0 h = {v e L 2 {VL) | for all K G % there is a K G R : v\ K = a K } (1) 

Pl,P13d piecewise linear continuous finite element (2d, 3d), the degrees of freedom are the 
vertices values. 

Pl h ={ve H\n) I \/K G %; v\ K G Pi } (2) 
We can see the description of the rest of the F.E.S. in the full documentation of FreeFem++. 

7 Boundary Condition 

We will see in this section how it's easy to define the boundary condition (B.C.) with 
FreeFem++, for more information about these B.C., we refer to the full documentation. 

7.1 Dirichlet B.C. 

To define Dirichlet B.C. on a border Yd C M like w|r d = /, we can proceed as on(gammad,u=f ) , 
where u is the unknown function in the problem. 

The meaning is for all degree of freedom % of the boundary referred by the label "gammad" , 
the diagonal term of the matrix an = tgv with the terrible giant value tgv (=10 30 by default) 
and the right hand side b[i] = "(II^r)[i]" x tgv, where the " (Hf l g)g[i] n is the boundary node 
value given by the interpolation of g. (We are solving here the linear system AX = B, where 
A = (an). , ■ . and B = (hi)- , „ ). 

\ l J 1 i=l..n\j=l..m \ l H=l..n / 

If u is a vector like u = (ul, u2) T and we have wl|r d = /I and «2|r d = /2, we can proceed 
as on (gammad, ul=f 1 ,u2=f 2) . 

7.2 Neumann B.C. 

du 

The Neumann B.C. on a border T n C R, like — |r = g, appear in the Weak formulation 

on n 

of the problem after integrating by parts, for example i = id' < ^ ) )r n = J 9 ' ®dx = 

intld(Th,gamman) (g*phi). 

7.3 Robin B.C. 

du 

The Robin B.C. on a border T r C M; like au + k— = b on F r where a = a(x,y) > 0, 

on 

y) > and b = b(x,y); also appear in the Weak formulation of the problem after 
integrating by parts, for example — = ( au ~ ^ < ^ ) )r = J au-<& dx — j b- $ dx = 
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intld(Th,gammar) (a*u*phi)-intld(Th,gammar) (b*phi). 

Important: it is not possible to write in the same integral the linear part and the bilinear 
part such as in intld(Th,gammar) (a*u*phi-b*phi) . 

7.4 Periodic B.C. 

In the case of Bi-Periodic B.C., they are achieved in the definition of the periodic F.E.S. 
such as : 

fespace Vh ( Th , PI , periodic = [ [1 , x] , [3 , x] , [2 , y] , [4 , y] ] ); 



(3) 



8 Solve the problem 

We present here different way to solve the Poisson equation : 
Find u : Q =]0, l[x]0, 1[ — ► R such that, for a given / e L 2 (fi): 

—Am = / in fl 
u = on dfl 

Then the basic variational formulation of (3) is : 
Find u G H£(Q), such that for all v G H^(Q), 

a(u, v) = l(v) (4) 

where 



a(u, v) — I V« • Vf dxdy and l(v) — / / • v dxdy 
Jn Jq 

To discretize (4), let Th denote a regular, quasi uniform triangulation of Q with triangles of 
maximum size h < 1, let Vh = {vh G C°(Q);v h \ T G Pi (T),VT G 7/,;^ = on dfl} denote a 
finite-dimensional subspace of Hq(Q) where Pi is the set of polynomials of R of degrees < 1. 
Thus the discretize weak formulation of (4) is : 



Find u h G V h : / Vu h ■ Vv h dxdy - I f ■ v h dxdy = Vv h G V h . 
Jn Jn 



(5) 



8.1 solve 

The first method to solve (5) is to declare and solve the problem at the same time by using 
the keyword solve such as : 

solve poisson (uh , vh , init=i , solver=LU) = // Solve Poisson Equation 

int2d(Th)( Grad (uh) ' *Grad (vh) ) // bilinear form 

-int2d (Th) (f *vh) // linear form 

+on (1 , 2 ,3 ,4 , uh=0) ; // Dirichlet B.C. 
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The solver used here is Gauss' LU factorization and when init 7^ the LU decomposition is 
reused so it is much faster after the first iteration. Note that if the mesh changes the matrix is 
reconstructed too. 

The default solver is sparsesolver ( it is equal to UMFPACK if not other sparce solver is 
defined) or is set to LU if no direct sparse solver is available. The storage mode of the matrix 
of the underlying linear system depends on the type of solver chosen; for LU the matrix is 
sky-line non symmetric, for Crout the matrix is sky-line symmetric, for Cholesky the matrix 
is sky-line symmetric positive definite, for CG the matrix is sparse symmetric positive, and for 
GMRES, sparsesolver or UMFPACK the matrix is just sparse. 



8.2 problem 

The second method to solve (5) is to declare the problem by using the keyword problem, 
and then solve it later by just call his name, such as : 

problem poisson (uh , vh , init=i , solver=LU) =// Definition of the 
^problem 

int2d(Th)( Grad (uh) ' *Grad ( vh) ) // bilinear form 

-int2d (Th) (f *vh) // linear form 

+on (1 , 2 ,3 ,4 , uh=0) ; // Dirichlet B.C. 

Poisson; // Solve Poisson Equation 

Note that, this technique is used when we have a time depend problem. 



8.3 varf 

In FreeFem++, it is possible to define variational forms, and use them to build matrices and 
vectors and store them to speed-up the script. 



The system (5) is equivalent to : 

Find u h e V h : a(u h , v h ) = l(v h ) Wv h e V h . (6) 

Here, 

M-1 

Uh(x,y) = ^2u hi (f)i(x,y) (7) 

i=0 

where 0, = Vhi,i = 0, M — 1 are the basis functions of Vh, M = Vh.ndof is the number of 
degree of freedom (i.e. the dimension of the space Vh) and um is the value of Uh on each degree 
of freedom (i.e. Uh% =uh[] [i]=£7). 



Thus, using (7), we can rewrite (6) such as : 

M-1 

A v u hj - Fi = 0, i = 0, • • • , M - 1; (8) 

3=0 

where 

: dxdy 



Aij = / V0j V0i dxdy and F = 

Jn Jn 
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The matrix A = (Aij) is called stiffness matrix. 

We deduce from the above notation that (8) is equivalent to 

A- U = F U = A' 1 ■ F (9) 

which can be solve in FreeFem++ as : 
int m=10,n=10; 

mesh Th= square (m,n ,[x,y]); 
f espace Vh ( Th , PI ) ; 
Vh uh , vh ; 

macro Grad (u) [dx (u) , dy (u) ] // in 2D 
func f = 1 ; 

varf a(uh ,vh) = int2d (Th) ( Grad (uh) '* Grad (vh) ) // bilinear 
f orm 

+on (1 ,2 ,3 ,4 , uh =0); // Dirichlet B.C. 
matrix A=a(Vh ,Vh); // build the matrix 

varf 1( unused ,vh) = int2d (Th)(f*vh); // linear form 
Vh F; F [] = 1(0, Vh) ; // build the right hand side vector 
set (A, solver = sparsesolver ); 
uh [] = A" -1*F [] ; 
plot (uh) ; 

And in 3D : 

load "msh3" 
load " medit " 
int m=10,n=10; 

mesh Th2= square (m,n , [x,y]); 

mesh3 Th=buildlayers (Th2 , 10 , zbound= [0 , 1] ) ; 

f espace Vh ( Th , P13d ) ; 

Vh uh , vh ; 

macro Grad (u) [dx (u) , dy (u) , dz (u) ] // in 2D 
func f=l; 

varf a(uh ,vh) = int3d (Th) ( Grad (uh) '* Grad (vh) ) // bilinear 
f orm 

+on (0, 1 ,2 ,3 ,4 ,5 , uh =0); // Dirichlet B.C. 
matrix A=a(Vh ,Vh); // build the matrix 

varf 1( unused ,vh) = int3d (Th)(f*vh); // linear form 
Vh F; F [] = 1(0, Vh) ; // build the right hand side vector 
set (A, solver = sparsesolver ); 
uh [] = A~ -1*F [] ; 
medit ( " sol " , Th , uh) ; 
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9 Learning by examples 



9.1 Rate of convergence for the Poisson equation 

At the beginning, we prove that the rate of convergence in space for the Poisson equation 
code with Pi finite element is of order 2. 

In this example, we took zero Dirichlet homogenous B.C. on the whole boundary and we have 
considered the following exact solution : 

u ex = sin(7rx) • sin(7n/) 

Then, we compute the corresponding right hand side f(x,y) in order to obtain the L 2 norm of 
the error between the exact solution and the numerical one (cf. Table 1) 

E(u, h n ) = \u h (h n ) - u ex (h n )\ L 2, Mi = 1, nref , h = 5x = 1/N; 
and then the rate of convergence in space 

= lo S (E(u,k JIE^)) = i ^ 

log (hn-x/hn) 



We give here a method to compute the right hand side using Maple 



6 



> u == una/pply (sm(pi-x) -sin(pi-j), x, y); 

u:= (x,y)^>sin(n x) sin(7ij) (1) 

> / : = -diff (\\{x,y),x,x) -diff (u(x, y), y, y); 

2 

/:= 2 sin(7i x ) 71 sin(Ttj) (2) 

> with ( CodeGeneration ) ; 

[ C, Fortran, IntermediateCode, Java, LanguageDefinition, Matlab, Names, Save, Translate, (3) 
VisualBasic\ 

> Matlab ( %%, resultname = " f ' ) ; 

f = 0.2el * sin (pi * x) * pi A 2 * sin (pi * y); 



We can copy and paste the result of f(x,y) in the FreeFem++ code. 
We present here the script to compute the rate of convergence in space of the code solving the 
Poisson equation : 

int nref=4; 

real[int] L2err or ( nref ) ; // initialize the L2 error array 

for (int n=0 ; n<nref ; n++) { 

int N=2~(n+4); // space discretization 

mesh Th= square(N,N); // mesh generation of a square 
fespace Vh(Th,Pl); // space of PI Finite Elements 
Vh uh , vh ; // uh and vh belongs to Vh 

6. http : //www.maplesof t . com/ 
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macro Grad (u) [dx (u) , dy (u) ] // 

Vh uex= sin (pi *x) * s in (pi *y ) ; // exact solution 
Vh f =0 . 2 e 1 * sin (pi *x) *pi "2* sin (pi *y ) ; // corresponding RHS 
varf a(uh,vh) = int2d(Th)( Grad (uh) ' *Grad ( vh) ) // bilinear 
orm 

+on (1 , 2 ,3 ,4 , uh=0) ; // Dirichlet B.C. 
matrix A=a(Vh,Vh); // build the matrix 

varf 1 (unused ,vh) = int2d (Th) (f *vh) ; // linear form 
Vh F; F [] = 1(0, Vh); // build the right hand side vector 

set (A , solver=sparsesolver ) ; 
uh [] = A~-1*F [] ; 

L2error[n]= sqrt ( int2d (Th) ( (uh-uex) "2) ) ; 

} 

for(int n=0 ; n<nref ; n++) 

cout << "L2error " << n << " = "<< L2error [n] <<endl; 
for(int n=l ; n<nref ; n++) 

cout <<" convergence rate = "<< log ( L2err or [n - 1] / L2err or [n] ) / 
**log(2.) <<endl; 



N n 


E(u, h n ) 


r(u, h n ) 


16 


0.0047854 




32 


0.00120952 


1.9842 


64 


0.000303212 


1.99604 


128 


7.58552e-05 


1.99901 



Table 1: L 2 norm of the error and the rate of convergence. 



9.2 Poisson equation over the Fila's face 

We present here a method to build a mesh from a photo using Photoshop® and a script in 
FreeFem++ made by Frederic Hecht. 

We choose here to apply this method on the Fila's face. To this end, we start by the photo in 
Figure 11, and using Photoshop®, we can remove the region that we wanted out of the domain 
such as in Figure 12, then fill in one color your domain and use some filter in Photoshop® in 
order to smooth the boundary as in Figure 13. Then convert your jpg photo to a pgm photo 
which can be read by FreeFem++ by using in a terminal window : 

convert fila.jpg fila.pgm 
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Figure 12: Using Photoshop. 
Figure 11: Initial photo. Figure 13: Last photo. 



Finally, using the following script : 



set to an rect 



array 



load "ppm2rnm" 
load " isoline " 
string f ila= " f i la . pgm " ; 
real [int , int] Curves (3,1); 
int [int ] be ( 1 ) ; 
int nc ; 

{ // build the curve file xy . txt . . . 
real [ int , int ] ffl(fila); // read image and 
// remark (0,0) is the upper, left corner, 
int nx = ffl.n, ny=ffl.m; 

// build a cartesain mesh such that the origne is qt the right 
^place . 

mesh Th = squar e (nx -1 , ny - 1 , [(nx-l)*(x) ,(ny-l)*(l-y)]) ; 
// warning the numbering is of the vertices (x,y) is 
// given by $ i = x/nx + nx* y/ny $ 
f espace Vh(Th , PI) ; 
Vh f 1 ; 

// transforme array in finite element function. 
= f 1 [] . max ; 
= f 1 [] . min ; 
( vmin+vmax ) /2 ; 



f 1 [] =f f 1 ; 

real vmax 
real vmin 
real vm = 
verbosity =3 ; 
/* 

Usage of isoline 
the named parameter : 
iso=0.25 // value of 
close=l, // to force 
beginend=be , // begin 
smoothing= . 01 , // nb 
where size is the size 



ISO 

to have closing curve . . . 
and end of curve 
of smoothing process = size 
of the curve . . . 



ratio * 0.01 



21 



ratio=0 . 5 

f ile = " f ilename " 

ouptut : 

xx, yy the array of point of the iso value 
a closed curve number n is 

in fortran notation the point of the curve are: 

(xx [i] , yy [i] , i = iO , il ) 

with : iO = be [2*n] , il = be [2*n + l] ; 

*/ 

nc=isoline (Th , fl , iso=vm , close=0 , Curves , beginend=be , smoothing 

■*■= . 005 , ratio=0 . 1) ; 
verbosity=l; 
} 

int ic0 = be(0), icl = be(l)-l; 

plot ( [Curves (0, icO : icl) , Curves (1 , icO : icl)] , wait = l) ; 
// end smoothing the curve .... 
macro GG(i) 
border G#i(t=0,l) 
{ 

P = Curve (Curves , be ( i *2) , be (i *2+l) -1 , t ) ; 
label = i + l ; 

} 

real lg#i=Curves (2 , be ( i *2+ 1 ) - 1) ; // 

GG(0) GG(1) GG(2) GG(3) GG(4) // number of closing curve 

real hh= -10; 

cout << " . . " <<endl ; 

func bord = GO ( IgO/hh) +G1 ( lgl /hh) +G2 ( lg2/hh) +G3 ( lg3/hh) +G4 ( lg4/hh 
*0 ; 

plot (bord , wait =1); 
mesh Th=buildmesh (bord) ; 
cout <<"... " <<endl ; 
plot (Th , wait = 1) ; 

Th = adaptmesh (Th , 5 . , IsMetric=l,nbvx=le6) ; 
plot (Th , wait=l) ; 
savemesh(Th, "fila.msh") ; 

we can create the mesh of our domain (cf. Figure 14), and then read this mesh in order to solve 
the Poisson equation on this domain (cf. Figure 15) 

mesh Th ( " f ila . msh " ) ; 
plot (Th) ; 

f espace Vh(Th , PI) ; 
Vh uh , vh ; 
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f unc f = 1 . ; 

macro Grad (u) [dx (u) , dy (u) ] // 

solve Poisson (uh , vh) = int2d (Th) (Grad (uh) ' *Grad (vh) ) - int2d(Th)( 

f*vh) + on (1 , 2 , 3 ,4 , 5 , uh = 0) ; 
plot(uh,dim=2,fill=true ,value=true) ; 




Figure 14: Mesh of the Fila's face. Figure 15: Solution on the Fila's face. 



9.3 Rate of convergence for an Elliptic non linear equation 

Let Q = B(0,R) C M 2 , it is proposed to solve numerically the problem which consist to find 
u(x, y) such that 

( -Au(x,y)+u 3 = f for all (x,y) G Q C IR 2 , , . 

\ u(x,y) = for all (x,y) on Oil. 

9.3.1 Space discretization 

Let Th be the triangulation of Q and 

V h = {v h E C°(n);v h \ T E Pi (T),VT G T h ,v h = on dil}. 

For simplicity, we denote by V(u) = u 2 , then the approximation of the variational formulation 
will be : 

Find Uh G Vh such that Wvh G Vh we have : 

- (Au h ; v h ) + (V(uh) ■ u h ] v h ) = (/; v h ) 

thus 

(Vu h ; Vv h ) + (V(u h ) ■ u h ; v h ) = (/; v h ) (11) 

In order to solve numerically the non linear term in (13), we will use a semi-implicit scheme 
such as : 

< V< +1 ; Vv h ) + (V«) • < +1 ; v h ) = (/; v h ) , 
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and then we solve our problem by the fixed point method in this way : 

Set u° h = u = 

Set VK) = V(tto) 

Set err = 1. 

while (err > le — 10) 

Solve (Vu p h +l ; Vv h ) = - (V«) -u p h ;v h ) + (f;v h ) 
Compute err = \\v% — 
set V«) = V« +1 ) 
set < = < +1 
p = p + 1; 

End while 

In order to test the convergence of this method we will study the rate of convergence in space 
(cf. Table 3) of the system (12) with R = 1 and the exact solution : 

u ex = sin(x 2 + y 2 - 1). 

Then, we compute the corresponding right hand side f(x,y) using Maple such as : 



> u == unapply (sin( x 1 + y 1 — l),x,y); 

u := (x, y) — >-sin(;r 2 + y 2 — 1 ) (1) 

> / : = -dif{u{x,y),x,x)-ffif (u{x,y),y,y) +u{x,y) 3 ; 

/:= 4 sin(,r 2 + y 2 — 1 ) x 2 — 4 cos(x 2 + y 2 — 1 ) +4 sin(x 2 — 1 ) y 2 + sin(x 2 + y 2 (2) 

-i) 3 

> simplify {%); 

4 sin(x 2 + y 2 — \ ) x 2 — 4 cos(x 2 + y — l) + 4 sm(x 2 + y 2 — l) y 2 + sin(x 2 + y 2 (3) 

— 1 ) — sin(x 2 + y 2 — l) cos ( x 2 + y 2 — 1 ) 

> with ( CodeGeneration ) ; 

[ £7, Fortran, IntermediateCode, Java, LanguageDefinition, Matlab, Names, Save, Translate, (4) 
VisualBasic\ 

> Matlab ( %%, resultname = "f ' ) ; 

f = 0.4el * sin((x A 2 + y A 2 - 1 ) ) * (x A 2 ) - 0.4el * cos((x A 
2 + y A 2 - 1)) + 0.4el * sin( (x A 2 + y A 2 - 1 ) ) * (y A 2 ) + sin 
((x A 2+y A 2-l))- sin((x A 2+y A 2-l)) * cos((x A 2 + y 
A 2 - 1 ) ) A 2 ; 



We present here the corresponding script to compute the rate of convergence in space of the 
code solving the Elliptic non linear equation (12): 

verbosity =0 . ; 
int nraff=7; 

real[int] L2err or ( nraf f ) ; // initialize the L2 error array 
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// space discretization 



for (int n=0 ; n<nraf f ; n++) { 
int N=2~ (n+4) ; 
real R=l.; // radius 

border C(t=0. ,2.*pi){x=R*cos(t) ;y=R*sin(t) ; label = 1 ; } ; 
mesh Th=buildmesh (C (N) ) ; 
f espace Vh(Th , PI) ; 

vh ; 
~ 2 
" 2 + 
- 1) ) 

((x " 



Vh 
Vh 
Vh 



2 

y 



y 
+ 

2 



D) ; 

" 2 - 
0.4el 

+ y " 

2 



uh , uhO=0, V=uhO 
uex=sin ( (x ~ 2 + 
f=0.4el * sin((x 
*-((x ~ 2 + y ~ 2 
(y " 2) + sin 

- 1) ) * cos((x * 2 + y * 2 - 1)) 
macro Grad (u) [dx (u) , dy (u) ] // 
problem ELLNL (uh , vh ) = 

int 2d (Th) (Grad (uh) ' *Grad (vh) ) 

+ int2d(Th) ( uh*V*vh ) 

- int2d (Th) ( f *vh ) 

+ on (1 , uh=0) ; 
real err=l . ; // for the convergence 
while (err >= le-lOH 
ELLNL ; 

err = sqrt ( int 2d (Th) ( (uh-uhO) "2) ) ; 
V=uh~2; // actualization 

uhO=uh ; 



1) ) * (x 
* sin ( (x 
2 - 1) ) ■ 



2) - 0.4el * cos 
2 + y ~ 2 - 1)) 
sin((x " 2 + y 



2; 



// bilinear term 

// non linear term 
// right hand side 
// Dirichlet B.C. 



L2error[n]= sqrt ( int 2d (Th) ( (uh-uex) "2) ) ; 

} 

for(int n=0 ; n<nraf f ; n++) 

cout << "L2error " << n << " = "<< L2error [n] <<endl; 
for(int n=l ; n<nraf f ; n++) 

cout <<" convergence rate = "<< log ( L2error [n - 1] / L2err or [n] ) / 
*»log(2.) <<endl; 



N n 


E(u, h n ) 


r(u, h n ) 


16 


0.015689 




32 


0.0042401 


1.88758 


64 


0.00117866 


1.84695 


128 


0.00032964 


1.83819 


256 


8.48012e-05 


1.95873 


512 


1.9631e-05 


2.11095 


1024 


4.88914e-06 


2.00548 



Table 2: L 2 norm of the error and the rate of convergence. 
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9.4 Rate of convergence for an Elliptic non linear equation with big 
Dirichlet B.C. 



(12) 



Let Q = B(0,R) C IR 2 , it is proposed to solve numerically the problem which consist to find 
u(x, y) such that 

Au(x, y) = V(u) ■ u for all (x, y) G Q C IR 2 , 
u(x,y) = DBC — >■ +00 for all (x, y) on dQ. 

9.4.1 Space discretization 

Let Th be the triangulation of Q and 

\4 = {v h G C°{Cl);v h \ T G Px(T),VT G %,v h =pon 8Q},p — >• +00. 

Then the approximation of the variational formulation will be : 
Find Uh G Vft such that Vf /, G 14 we have : 

(Au h ;v h ) = (V(u h ) ■ u h ;v h ) 

thus 

- (V« h ; Vv ft ) = (V(« h ) • u h ) (13) 

In order to solve numerically the non linear term in (13), we will use a semi-implicit scheme 
such as : 

-<V< +1 ;V^> = <V«)-< +1 ;^>, 
and then we solve our problem by the fixed point method in this way : 

Set u ° h = u = DBC,p = 

Set V«) = V(«o) 

Set err = 1. 

while (err > le — 10) 

Solve - (V< +1 ; Vu h > = (V«) • <; u h ) 
Compute err = \\v% — u 1 ^ 1 \\l2 
set V«) = V« +1 ) 
p = p + 1; 



End while 



In order to test the convergence of this method we will study the rate of convergence in space 
(cf. Table 3) for an application of (12), where R = 1, V(u) = u, DBC = or DBC = 50. 
The system to be solved is then 

Au(x, y)-u 2 = f for all (x, y) G Q = B(0, 1) C 1R 2 , (14) 
u(x,y) = DBC for all (x, y) on dQ. (15) 

In this case, we will the following exact solution : 

u ex = DBC + sin(x 2 + y 2 - 1). 

Then, we compute the corresponding right hand side f(x,y) using Maple such as : 

We present here the corresponding script to compute the rate of convergence in space of the 
code solving the Elliptic non linear equation (14): 
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> u == unapply ( DBC + sin ( x 2 + y 2 — l),x,y); 

u\= {x,y)^DBC + sin(x 2 + y 2 - l) (1) 

> /'■= diff (u(x,y), x, x) + diff (u(x, y ), y, y ) - u(x, y) 2 ; 

/:= -4sin(x 2 - \ ) x 2 +4cos(x 2 +/ - l) -4sin(x 2 +/ - l) / - (2) 

+ sm{x 2 + y 2 — \ )) 

> with ( CodeGeneration ) ; 

[ £7, Fortran, IntermediateCode, Java, LanguageDefinition, Matlab, Names, Save, Translate, (3) 
VisualBasic\ 

> Matlab ( %%, resultname = "f ) ; 

f = -0.4el * sin((x " 2 + y A 2 - 1 ) ) * (x A 2 ) + 0.4el * cos((x " 
2 + y A 2 - 1)) - 0.4el * sin( (x A 2 + y " 2 - 1 ) ) * (y A 2 ) - 
_(DBC + sin((x "2+y"2-l)))"2; 



initialize the L2 error array 



// space discretization 



int nref=7; 

real[int] L2err or ( nref ) ; // 
for (int n=0 ; n<nref ; n++) { 

int N=2~ (n+4) ; 

real R=l.; // radius 

border C(t=0. ,2.*pi){ x = R*cos(t) ;y = R*sin(t) ; label = 1 ; } ; 
mesh Th=buildmesh (C (N) ) ; 
f espace Vh(Th ,P1) ; 



real DBC=0 . ; 

Vh uh, uhO=DBC, V=uhO , vh 
Vh uex=DBC+sin ( (x " 2 + y 
Vh f=-0.4el * sin((x " 2 
*»cos((x " 2 + y " 2 - 
*•!)) * (y ~ 2) - (DBC 



" 2 - 1) ) ; 
+ y ~ 2 - 1)) 
1) ) - 0.4el * 
+ sin((x ~ 2 ■ 



* (x " 
s in ( ( x 



2) + 0.4el * 
~ 2 + y ~ 2 
- 1))) ~ 2; 



macro Grad (u) [dx (u) , dy (u) ] // 
problem ELLNL(uh,vh) = 

- int2d(Th) (Grad(uh) '*Grad(vh)) 

- int2d(Th) ( uh*V*vh ) 

- int2d (Th) ( f *vh ) 
+ on (1 , uh = DBC) ; 

real err=l.; // for the convergence 
while (err >= le-10){ 
ELLNL ; 

err = sqrt ( int 2d (Th) ( (uh-V) "2) ) ; 
V=uh; // actualization 

} 

L2error [n] = sqrt ( int 2d (Th) ( (uh-uex) "2) ) 



// bilinear term 

// non linear term 

// right hand side 

// Dirichlet B.C. 



for(int n=0 ; n<nref ; n++) 
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cout << "L2error " << n << " = "<< L2error [n] <<endl ; 
for(int n=l ; n<nref ; n++) 

cout <<" convergence rate = "<< log ( L2err or [n - 1] / L2err or [n] ) / 
**log(2.) <<endl; 



N n 


E(u,h n ), DBC=0 


r(u, h n ), DBC=0 


E(u, h n ), DBC=50 


r(u,h n ), DBC=50 


16 


0.0159388 




0.00610357 




32 


0.00455562 


1.80683 


0.00244016 


1.32268 


64 


0.00118025 


1.94855 


0.000767999 


1.6678 


128 


0.000335335 


1.81542 


0.000210938 


1.86429 


256 


8.6533e-05 


1.95428 


5.67798e-05 


1.89337 


512 


1.9715e-05 


2.13395 


1.40771e-05 


2.01203 


1024 


4.90847e-06 


2.00595 


3.56437e-06 


1.98163 



Table 3: L 2 norm of the error and the rate of convergence. 



9.5 Rate of convergence for the Heat equation 

Let Q =]0, 1[ 2 , we want to solve the Heat equation : 
du 

— —fi-Au = f(x,y,t) for all (x,y) GOcR 2 ,^? 



(16) 



u(x,y,0) = u (x,y), 

u = on dfl. 

9.5.1 Space discretization 

Let Th be the triangulation of Q and 

V h = {v h G C°(tiy,v h \ T G Px(T),VT G %,v h = on dtt}. 

Then the approximation of the variational formulation will be : 
Find Uh G Vh such that Wft G Vh we have : 

v h ^j - (ji ■ Au h ; v h ) = (/; v h ) 
Thus 

-^r;v h )+v (VM h ; Vv h ) = (/; v h ) (17) 

9.5.2 Time discretization 

We will use here a ^-scheme to discretize the Heat equation (17) as : 

n+l _ „,n 

At ' 



u h - u h 



Vh\+n-6(Vu n h +1 ;Vv h )+n-(l-6) (V<; Vv h ) = 9 ■ (f n+1 ;v h ) + (1-9)- (f n ; v h ) . 
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Therefore 

,n+l 



U, 



; v h \+fi-9 <V< +1 ; Vv h ) = uA-^l-fl) (V<; V^)+0-(.T +1 ; %)+(l-£)-(r ; Uft > 

(18) 



At 

To resolve (18) with FreeFem++, we will write it as a linear system of the form : 



M-l 

AX = B^J2 Ai • Xj = Bi for i = 0; . . . ; M - 1; 

where, M is the degree of freedom, the matrix Aij and the arrays Xj and Bj are defining as : 

tgv = 10 30 if z G dQ and j = % 

,n+l 



t^u = 10 30 if z G <9fi 

^ = ^ / ~ ^ ' ( X ~ ^) • V < V ^* + ( ■ + ( x - 0) ■ f n ) Vidxdy otherwise 

We note that the ^-scheme is stable under the CFL condition (when 9 G [0, 1/2 [) : 

At At 1 



Ax ^Ay ~ 2- (1-26) 

In our test, we will consider that Ax = Ay and that CFL g]0, 1], then when 9 G [0, l/2[, the 
^-scheme is stable under this condition : 

CFL ■ (Ax) 2 



At < 



A -/I- (1 -29)' 



and for 9 G [1/2, 1], the ^-scheme is always stable. 
We note also that due to the consistency error : 

e\ 3 < cAt \29 - 1| + O (Ax 2 ) + O (Ay 2 ) + O (At 2 ) , 

the ^-scheme is consistent of order 1 in time and 2 in space when 9 = (with the CFL condition) 
and when 9 = 1 (with At = (Ax) 2 ) and the ^-scheme is consistent of order 2 in time and in 
space when 9 = 1/2 (with At = Ax) (cf. Table 4). 

Remark. When we use finite element, mass lumping is usual with explicit time-integration 
schemes (as when 9 G [0, l/2[). It yields an easy-to-invert mass matrix at each time step, 
while improving the CFL condition [Hug87]. In FreeFem++, mass lumping are defined as 
int2d (Th , qf t =qf lpTlump) . 

In order to test numerically the rate of convergence in space and in time of the ^-scheme, we 
will consider the following exact solution : 

u ex = sin(7rx) ■ sin(7n/)e sm(i) . 

Then, we compute the corresponding right hand side f(x,y) using Maple such as : 
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> u == unapply (sin(pi-x ) -sin(pi-j) -exp(sin(/) ), x,y, /); 

u '.= (x, y, /) — >-sin(7i x ) sin(7l>') e sm ^' (1) 

> / : = diff {u{x, y, t), /) - vmx-diff (u{x,y, t), x, x) - vmx-diff (u{x, y, t), y, y ); 

_/*:= sin(7i sin(7i cos(/) + 2 |j, sin(7ix) 71 sin(xc^) e sm ^' (2) 

> simplify {%); 

sin(jcx) sm{ny) e sm(0 (cos(/) +2\in 2 ) (3) 

> with ( CodeGeneration ) ; 

[ C, Fortran, IntermediateCode, Java, LanguageDefinition, Matlab, Names, Save, Translate, (4) 
VisualBasic\ 

> Matlab ( %%, resultname = "f ) ; 

f = sin (pi * x) * sin (pi * y) * exp(sin(t)) * (cos(t) + 0.2el * mu 
_* Pi A 2); 



we remind that the rate of convergence in time is 

r(M«„,M = M^plV|(^M , V „ = 1 nref 
log (dt n _i/dt n ) 

macro Grad (u) [dx (u) , dy (u) ] // 

macro uex(t) ( sin (pi *x) * sin (pi *y ) * exp ( sin ( t ) ) ) // 

macro f(t) ( sin(pi * x) * sin(pi * y) * exp(sin(t)) * (cos(t) + 
0.2el * mu * pi ~ 2) ) // 



real t, dt , h,T=.l, mu=l., CFL=1., theta=0.; 
int nref=4; 

real[int] L2err or ( nref ) ; // initialize the L2 error array 

real [int] Dx(nref); // initialize the Space discretization array 

real [int] DT(nref); // initialize the Time discretization array 



for (int n=0 ; n<nref ; n++) { 
int N=2~ (n+4) ; 
t = ; 
h=l ./N; 
Dx [n] =h ; 
if (theta<.5) 

dt=CFL*h~2/4 . / ( 1 . -2 . * theta) /mu ; 
else if (theta==.5) 

dt = h ; 

else 

dt=h~2 ; 
DT [n] =dt ; 

mesh Th=square (N , N) ; 
f espace Vh(Th , PI) ; 
Vh u , uO , B ; 
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varf a(u,v) = int2d (Th , qf t = qf lpTlump ) (u* v/dt + Grad (u) ' * Grad ( 

**v) *theta*mu) + on ( 1 , 2 , 3 , 4 , u = 0) ; 
matrix A = a(Vh,Vh); 

varf b(u,v) = int2d (Th , qf t =qf lpTlump ) (uO* v/dt - Grad(uO)'* 
^•Grad (v) * (1 . -theta) *mu ) + int2d (Th , qft = qf lpTlump ) ( (f(t + 
**dt) *theta + f (t) * (1 . -theta) ) *v ) + on ( 1 , 2 , 3 , 4 , u = 0) ; 

u=uex (t ) ; 

for (t=0 ; t<=T ; t+=dt ) { 
uO = u ; 

B [] = b(0 , Vh) ; 

set (A , solver = spar se solver ) ; 

u [] = A"-1*B [] ; 

} 

L2error [n] = sqrt ( int2d (Th) ( abs (u-uex ( t ) ) " 2) ) ; 

} 

for(int n=0 ; n<nref ; n++) 

cout << "L2error " << n << " = "<< L2error [n] <<endl; 
for(int n=l ; n<nref ; n++) { 

cout <<"Space convergence rate = "<< log ( L2error [n-1] /L2error 

*-[n] ) /log(Dx [n-1] /Dx [n] ) <<endl; 
cout <<"Time convergence rate = "<< log ( L2error [n- 1] / L2error [ 
*»n] ) /log (DT [n-1] /DT [n] ) <<endl; 

} 



N n 


16 


32 


64 


128 


E( u ,h n ),e = o 


0.00325837 


0.000815303 


0.000203872 


5.09709e-05 


r(u, h n ), 9 = 




1.99874 


1.99967 


1.99992 


r(u,dt n ,h n ), 9 = 




0.99937 


0.999834 


0.99996 


E(u,h n ),9 = l/2 


0.00325537 


0.000819141 


0.000203817 


5.08854e-05 


r(u,h n ), 9 = 1/2 




1.99064 


2.00684 


2.00195 


r(u,dt n ,h n ), 9 = 1/2 




1.99064 


2.00684 


2.00195 


E{u,h n ), 9 = 1 


0.00323818 


0.000807805 


0.000201833 


5.04512e-05 


r(u, h n ), 9 = 1 




2.0031 


2.00084 


2.0002 


r(u,dt n ,h n ), 9 = 1 




1.00155 


1.00042 


1.0001 



Table 4: L 2 norm of the error and the rate of convergence in space and in time for different 9. 



10 Conclusion 

We presented here a basic introduction to FreeFem++ for the beginner with FreeFem++. For 
more information, go to the following link http: //www. f reef em.org/ff++. 
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